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We reveal that the transition in the saturated flux for aeolian saltation is generically discontinuous 
by explicitly simulating particle motion in turbulent flow. This is the first time that a jump in the 
saturated flux has been observed. The discontinuity is followed by a coexistence interval with two 
metastable solutions. The modification of the wind profile due to momentum exchange exhibits a 
maximum at high shear strength. 



Aeolian saltation is one of the main actors of molding 
Earth's landscape. It not only has impact on the evo- 
lution of the agricultural areas but also on the change 
of river courses and the motion of sand dunes. Salta- 
tion consists in the transport of sand or gravel by air or 
water and occurs when sand grains are lifted and accel- 
erated by the fluid, hopping over the surface and ejecting 
other particles [l|, 0] ■ The proper measurement of salta- 
tion close to the impact threshold through experiments 
is strongly limited by technical difiiculties. Because of 
temporal fluctuations, one obtains unreliable data and 
large error bars. The computational simulation of ae- 
olian saltation, however, can monitor every mechanical 
interaction between particles giving local insight about 
the particle splash and, in particular, about the system 
close to the impact threshold. Here, we report for the 
first time the existence of a discontinuity in the flux at 
the onset of saltation using discrete elements. 

Saltation was first described by Bagnold's seminal 
work discussing the particle splash and proposing a cubic 
relation between wind shear velocity and saturated mass 
flux Q-Q. Greeley et al. Q confirmed experimentally 
Bagnold's result and Ungar and Haff [3] complemented it 
in their numerical model with a quadratic relation near 
the impact threshold. Anderson and Haff 0, Q stud- 
ied numerically the statistical properties of the particle 
splash relating the number and velocities of ejected par- 
ticles to the impacting ones. Splash entrainment was 
carefully studied in experiments |10l4l3j and used in non- 
steady saltation models [3, [l^ finding that the splash 
mechanism dominates the saturation of sand flux. Re- 
cently, Almeida et al [l^, [13] implemented the full feed- 
back with the fluid and flxed splash angle while Kok and 
Renno [Tsj used a splash entrainment function, both re- 
sulting in good agreement with experiments. These and 
other theoretical studies, e.g., |19l - [2l| . however, describe 
the sand bed as a rough wall instead of resolving it at 
the particle scale and consequently rely on an empirical 
splash function. 
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FIG. 1. [Color online] Snapshot of the simulation with an ini- 
tially logarithmic wind profile. Colors represent the particle 
velocity. 
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where j/q is the roughness of the bed, /iq is the bed height, 
K — 0.4 the von Karman constant, and the wind shear 
velocity. Below Hq, the wind velocity is zero. We adopted 
the widely used roughness law, yo = Dmean/30, mea- 
sured by Refs. [13, [Ijl for hydrodynamically rough pipe 
flow and conflrmed by Ref [5] for air flow, where D^ean 
is the mean diameter of the particles. A particle of di- 
ameter D is accelerated by the wind drag force given by 

a 



We present a discrete element simulation for aeolian 
saltation which does not require the use of a splash func- 
tion. Particles are subjected to gravity g and dragged by 
a height-dependent wind field. The unperturbed wind 
profile is [3] 



PaCdVrVr 
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where pa is the air density and Vr = v — u is the velocity 
difference between particle and air, with u,. = |vr|. The 
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drag coefficient Cd proposed by Cheng [241 is suited to 
model natural and irregularly shaped grains: 
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where /i ~ 1.8702 x 10^^kg/{m.s) is the dynamic viscos- 
ity and Re is the Reynolds number, which together with 
the Shields number 6 below are the pertinent dimension- 
less parameters of our problem: 
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where s — Ps/ Pw is the ratio between the grain and fluid 
density. 

The detailed integration of the wind profile considering 
the momentum exchange is explained in the Supplemen- 
tal Material [11]. 

An aerodynamic lift arises from shear in the fiow, 
which results in a pressure gradient normal to the shear 
in the direction of decreasing velocity. This aerodynamic 
lift can be approximately described by 
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where the lift coefficient C/ is proportional to Cd |26l |. 

The vertical motion of the particle is given by the com- 
petition between the gravity lift forces, and the re- 
bouncing of particles with the ground. 

Alternatively, instead of lift forces, we also did some 
simulations perturbating the system at rest by lifting up 
at every second a fraction c = 0.2 of the surface particles 
by a height Dmean with a probability c = 0.2 in order to 
restart the saltation. 

We consider a disordered particle bed with 500 spheri- 
cal particles initially at rest in a two dimensional system. 
The diameters are randomly chosen from a Gaussian dis- 



tribution around Dn 



of width 0.15L»„ 



We simu- 



lated systems with more particles to verify that the bot- 
tom wall effects are negligible as discussed in Ref [13, UH ■ 
Within the error bars, they displayed identical properties 
and therefore we did not need to consider larger systems. 
At i = 0, some particles are dropped from randomly cho- 
sen heights and when they reach the ground they collide 
with particles at rest inside the bed and thereby trigger 
saltation. The particle collisions are computed using the 
discrete elements method. 

Trajectories are obtained by integrating the equa- 
tions of motion according to the velocity-Stormer-Verlet 
scheme [IJ, using a spring dashpot potential with the 
spring constant k and a dissipative damping parameter 7. 
Further details about the technique are explained in the 
Supplemental Material [25|. The dissipative rate of the 
lower boundary was set high enough, 7^, = 0.8, to avoid 
that the shock wave generated by the impact reaches the 



system boundaries and be reflected to the bed surface, 
ejecting additional grains [ll|- 

Space is sliced in vertical rectangular domains of size 
(250 X 75)D,nean- The top is placed sufficiently high to 
mimic an open system. Periodic boundary conditions are 
imposed in wind direction and top and bottom bound- 
aries are reflective. 




FIG. 2. Flux series for two different Shields numbers near the 
onset. In (a), the system was perturbed at t = 9, 10, and lis 
to keep the flux diflferent from zero. The small peaks at t = 
9 and 10s show that some particles were lifted but did not 
provoke a splash. In (6) the system never settled down. 

We measure the particle flux through 
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where vf is the particle velocity in the x direction. The 
saturated flux is the average flux in the stationary state. 
Simulated systems with different number of particles dis- 
play similar flux. The dimensionless flux can be obtained 

by 



q = 



Ps^{s-l)gDf, 



(7) 



where ps is the grain density. 

We verified numerically that the dimensionless satu- 
rated fiux remains invariant under changes of the param- 
eters as long as the Reynolds number and Shields number 
are fixed. 

Simulations verify the existence of an onset velocity 
for sustained fiux 9c = 0.048 with strong temporal fluc- 
tuations as illustrated in Fig. O Below this value the 
flux will stop after some time. This happens because an 
impact may not necessarily eject other particles that can 
carry on saltation, and once the flux is zero it cannot 
restart unless a perturbation or lift force is introduced. 
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FIG. 3. [Color online] Saturated flux as function of the Shieds 
number, (a) Simulated data fitted by Qs = go + A{9 — 9c)0^^^ , 
qo = 0.01, A = 1.52, and 9c — 0.048 in comparison with 
experimental data ^2^.(6) Detailed view of the metastable re- 
gion of the discontinuous transition. The green curve with 
triangles and blue lines with squares are the results including 
perturbation with c = 0.2 and lift forces with Ci — 0.425C£j, 
respectively, (c) Average transient time to settle without 
perturbation in the metastable region. The dashed line fits 
the time distribution according to ts = to[9/{9c — 9)]^ with 
to = 0.524s, p = 1.5, 9c = 0.048. 



In Figl2lja) below the onset, some perturbations were 
introduced at t = 9, 10, and lis in order to trigger salta- 
tion, after it stopped. The perturbation at t = lis 
restarted the saltation and the flux returned to the pre- 
vious levels. For 9t < < 9c, the system displays this 
kind of metastable behavior with two possible solutions, 
either saltation or no motion, strongly dependent on the 
initial conditions and the triggering mechanism. 9t is the 
Shields number of the threshold, below which no particle 
is transported and 6c is a critical Shields number sepa- 
rating the metastable region from the normal regime of 
saltation. The previously defined perturbations or lift 
forces according to Eq. [5] are not sufficient to restart 
saltation for 9 < 9t, where 9t ~ 0.037 for perturba- 
tions with c = 0.2 and 9t — 0.036 for lift forces with 
Ci — 0.425Cd. The lower bound of the metastable region 
changes depending on the perturbation probability or the 
lift constant. They are no longer needed to maintain the 
saltation for 9 > 9c- The fact that perturbations are nec- 
essary to sustain saltation underlines the importance of 
the turbulent lift forces in the metastable region. 

Fig. Ela) presents the saturated flux for different 
Shields numbers in comparison with held and wind tun- 
nel experimental data from Iversen and Rasmussen [28j . 
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FIG. 4. [Color online] Momentum exchange in y direction 
becomes nonmonotonic as function of height as we increase 
the Shield number. In the bottom part, it increases with 
the wind velocity due to the high concentration of particles. 
The concentration shown in the inset decays logarithmically 
for lower velocities but has a pronounced inflection point for 
higher velocities. 



The calculated data are well fitted using = go + A{9 — 
9c)9^^^ with A = 1.52 and qo = 0.01, which is similar to 
the relation proposed in Ref.fl^. Figure IS^b) presents 
details of the discontinuous transition at 9c with a jump 
go in the saturated flux. The green triangles are the data 
including perturbations and the blue squares are the data 
with lift forces. Additional simulations show that 9t and 
the jump go depend on the lift coefficient Ci. 

Figure [3ljc) presents the average transient time <ts> 
it takes for the system to settle without perturbation or 
lift diverging when 9 ^ 9c- The dashed line fits the data 
according to = to[9/{9c - 9)Y with to = 0.052s, p = 
1.5. The transient times < tg > were averaged over 60 
simulations with different initial conditions. 

The discontinuity in the saturated flux is verified also 
at the same value 9c using the drag coefficient from Ref. 
[30| and is shifted in the velocity axis by the dissipation 
rate 7, i.e., to lower (higher) critical velocities for lower 
(higher) dissipation rates. Interestingly, we were able to 
find the same discontinuous transition and metastable 
region also using the code of Ref [iSj . 

Figure HI shows the momentum transfer to the particles 
as a function of height for different Shields numbers. Be- 
cause of the higher concentration of particles, the largest 
momentum transfer occurs close to the bed surface. Ad- 
ditionally, this region coincides with the one in which we 
observe small changes in the modified wind profile in Fig. 
1 in the Supplemental Material [25]. Momentum trans- 
fer for high Shields number exhibits a local maximum 
which gets more pronounced with increasing the Shields 
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FIG. 5. [Color online] Modified wind profiles intersect for 
6 > 9c a.t the Bagnold focus. This estimate corroborates 
measurements for a particle size Dmean ~ 0, HE, E, 



number. As shown in Fig. IH^b), high Shields numbers 
enhance the concentration at higher y where the wind 



profile grows stronger. This increases momentum trans- 
fer at this level. In Fig. IH^b), we see that in fact the 
concentration exhibits a saddle point at the same height 
at which the maximum in momentum transfer appears. 

Figure [5] shows modified velocity profiles for 9 > 9c 
crossingeach other at a focal point, called the "Bagnold 
focus" '5"] and approximately located 0.5cm above the 
sand bed, which is in qualitative agreement with mea- 
surements and theory 0, [H, [H, Hlj . 

In summary, we performed particle simulations with- 
out any assumption concerning particle trajectories or 
splash and revealed that the transition is discontinuous 
with a metastable state for aeolian saltation never re- 
ported before. In the metastable region, perturbations 
or lift forces are required to keep saltation going. It 
would be interesting to reanalyze or remake experiments 
of saltation very close to the onset to verify the predicted 
jump. The existence of a discontinuity at the threshold 
of turbulent particle transport can have far-reaching con- 
sequences in numerous applications ranging from dune 
mitigation or pneumatic transport to the understanding 
of Martian landscapes. 
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I. THE WIND PROFILE 

The description of the wind considers a logarithmic 
profile given by 



y-ho 

uiv) = — m . 



(1) 



But when particles are accelerated also the wind is 
decelerated, modifying the wind profile [T|. We obtain 
the grain stress profile Tg (y) by integrating the drag forces 
acting on all particles above height y, as 



oo 

,iy) = J f{y')dy' ^ ^ 
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where f{y') is the average force per unit volume and 
the drag force applied on the particles at a position 
above y, which depends on the wind velocity Uj at the 
same position, and A is the cross-sectional area parallel to 
the ground assuming a finite system thickness. The grain 
stress profile is then used to calculate the modified wind 
profile by rewriting Eq. ((ij as a differential equation. 



du 
dy 



Ky 



(3) 



where the modified wind strength Ur (y) depends accord- 
ing to Anderson and Half [1| on the grain stress (Eq. [2]) 
as 
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FIG. 1. [Color online] Wind profiles with and without momen- 
tum exchange. The modified profile has a weaker increase at 
highly concentrated regions but an asymptotic slope equal to 
the unmodified one. At y — ho, we have u = 0. 



Figure [T] shows an example of the wind profile for — 
1.4 m/s. As the wind strength decays with the number 
of saltating particles, the slope of the modified curve falls 
much below the original profile in the region nearby the 
particle bed. Because few or no particles are located 
in the higher regions, the slope of the modified curve 
asymptotically converges to the unmodified one. 
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The numerical solution is achieved iteratively. If no 
particle is accelerated Mt(j/) = w* and consequently, the 
unperturbed logarithmic profile is obtained. Otherwise, 
a grain stress arises reducing the shear velocity in the 
underneath areas. 

The wind profile is recalculated at each iteration step 
starting from the top of the bed for which basis position 
ho and shape change dynamically. Consequently, we need 
to recalculate also ho each time. Highly concentrated 
areas like the one close to the sand bed reduce strongly 
wind velocities. If the calculated velocity Vi in area yi is 
below O.lu*, this area always contains bed particles and 
the velocity is set to zero. This means that Hq is chosen to 
be the point where the calculated velocity exceeds O.lu*. 



II. PARTICE DYNAMICS 

Molecular dynamics is used to solve numerically the 
Newton's equations of motion. The differential equations 
are iteratively solved using time-discretization obtaining 
the new positions and velocities of the particles from the 
old positions, old velocities, and the corresponding forces. 
An efficient and, at the same time, stable approach for 
the time discretization of Newton's equations is the Ver- 
let algorithm which builds on the integration method of 
Stormer found in Ref Q. Let , vf and Ff be the 
vectors corresponding to the position, velocity and force 
of the particle i at time step n. Velocity and position of 
particle i at the next time step can be obtained by Q 
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where fc is a spring constant, is the mass of particle i, 
and Tij points from particle i to j and a dissipation due 
to the inelasticity of the collision, 



where St and rrii are the time step and the particle mass. 
Fjis the resultant force of all particles acting on particle 
i i 
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F.= E (7) 

When two particles i and j overlap (i.e. when their 
distance is smaller than the sum of their radia) two forces 
act on the particle i, an elastic restoration force, 

i^«=fcm4|r.,|-l/2K + d,)]^ (8) 



where Vjj = — is the relative velocity and 7 is phe- 
nomenological dissipation coefficient. In our work, we 
chose fc = 0.5 and 7 — 0.3 for particle collisions. The 
coefficient of restitution is given by the ratio between 
the velocities after and before the collision. We neglect 
Coulomb friction and rotation of particles. When a par- 
ticle collides with a wall the same forces act as if it would 
have encountered another particle of diameter do at the 
collision point 
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